# -*- coding: utf-8 -*-
"""
Created on Sun Jan 14 20:40:29 2024

@author: xiang
"""

import numpy as np
import matplotlib.pyplot as plt



def identity_plasma_init(self, flux_plasma, flux_limiter, init_factor=10, gapx=0):
    """
    find LCFS of plasma region with flux_plasma and flux_limiter. 
    judge plasma is limiter or divertor configuration with inductance on limiter

    input:
        option name         data_type    size     unit    meaning
        1      flux_plasma  array        (nr*nz,1)        flux of plasma induced by plasma and pfs
        1      flux_limiter array        (nlimiter,1)     flux of limiter induced by plasma and pfs
        0      isplot       int          1                is plot or not. 

    output:
        option  name            data_type   size    unit    meaning
        1      flux_plasma      array       (nr*nz,1)       flux of plasma induced by plasma and pfs
        1      flux_plasma_norm array       (nr*nz,1)       flux of plasma normalized
        1      index            array       (nr*nz,1)       index of plasma: 1 means inside LCFS; 0 means outside LCFS
        
    reference:
        [1]

    Developed by xiangguchn 2023/12/25
    """
            
    # plasma and limiter region
    x_plasma = self.eq['green']['grid']['plasma']['x']
    y_plasma = self.eq['green']['grid']['plasma']['y']
    
    x_limiter = self.eq['machine']['limiter']['x_limiter']
    y_limiter = self.eq['machine']['limiter']['y_limiter']
    # num_limiter = len(x_limiter)
    # num_init = len(x_init)
    
    if gapx==0:
        gapx = (x_plasma[1]-x_plasma[0])/2
    
    
    v_step = 1e-3 * init_factor
    ok_step = 5e-6 * init_factor
    
    my_delay = 0.5
    
    max_flux_plasma = max(flux_plasma)
    index_max_flux_plasma = np.argmax(flux_plasma)
    
    
    # plasma and limiter region
    
    # position higher than 1m should be cancelled
    y_top_limite = 1
    x_limiter_2 = x_limiter[abs(y_limiter)<y_top_limite]
    y_limiter_2 = y_limiter[abs(y_limiter)<y_top_limite]
    flux_limiter_2 = flux_limiter[abs(y_limiter)<y_top_limite]
    
    
    # find the max flux at limiter
    max_flux_limiter = max(flux_limiter_2)
    index_flux_limiter = np.argmax(flux_limiter_2)
    
    
    # find LCFS
    d_flux1 = max_flux_limiter
    d_flux0 = max_flux_plasma
    d_flux = d_flux1
    # plt.figure()
    # self.plot_flux(flux_plasma)
    while abs(d_flux0-d_flux1) > ok_step/10:
        d_flux = (d_flux1 + d_flux0)/2
        c = plt.contour(np.reshape(x_plasma, (129,129)), 
                        np.reshape(y_plasma, (129,129)), 
                        np.reshape(flux_plasma, (129,129)), 
                        levels=[d_flux])
        dp, cs = [], []
        for item in c.collections[0].get_paths():
            p = item.vertices
            dp.append((p[0,0]-p[-1,0])**2 + (p[0,1]-p[-1,1])**2)                
            cs.append(p)
            plt.plot(p[:,0], p[:,1], '.')
        
        d = min(dp)
        if d <  1e-5:
            d_flux0 = d_flux
        else:
            d_flux1 = d_flux
        
        index_d = np.argmin(dp)
        c_lcfs = cs[index_d]
        
    
    
    # find plasma inside lcfs
    rmin = min(c_lcfs[:, 0])
    rmax = max(c_lcfs[:, 0])
    zmin = min(c_lcfs[:, 1])
    zmax = max(c_lcfs[:, 1])
    self.fgs['xp'] = (rmin+rmax)/2
    
    index = flux_plasma*0
    index[flux_plasma>d_flux] = 1         
    index[x_plasma>rmax] = 0       
    index[x_plasma<rmin] = 0       
    index[y_plasma>zmax] = 0      
    index[y_plasma<zmin] = 0 
    
    # set plasma geo parameter
    self.fgs['Xp'] = (rmin+rmax)/2
    self.fgs['Yp'] = (zmin+zmax)/2
            
    flux_plasma_norm = (flux_plasma-max_flux_plasma)/(d_flux-max_flux_plasma)
    
    return flux_plasma, flux_plasma_norm, index
    # return flux_plasma_norm, index
